详解QT中的快速正弦、快速余弦函数的实现

  数值计算中，许多数学函数(如正弦、余弦、平方根等)可以通过泰勒级数展开的方法求得，然而，对于正弦函数与余弦函数，采用泰勒级数展开的方法收敛速度慢，计算速度慢。那么，有没有快速计算正弦函数与余弦函数并且可以估算误差限的算法呢？   QT中提供了快速正弦函数、快速余弦函数的算法，代码片段如下：

  QT中提供的快速正弦、快速余弦的算法，是一种典型的以空间换时间的查表法。这里以快速正弦函数为例进行讲解，对于快速余弦函数很容易类比推导得到。   由于 s i n ( x ) sin(x) sin(x)为周期为 2 π 2\pi 2π的函数，对于任意实数 x x x均可以变换到区间 [ 0 , 2 π ] [0,2\pi] [0,2π]，假设区间 [ 0 , 2 π ] [0,2\pi] [0,2π]的正弦表 s i n e T a b l e sineTable sineTable长度为 N N N。   求 s i n ( x ) sin(x) sin(x)时，实数 x x x所对应的最接近其值的正弦表索引(不一定在 0 至 N − 1 0 至 N-1 0至N−1之间): s i = r o u n d ( x N 2 π ) (1) s_i=round(\dfrac{xN}{2\pi})\tag 1 si​=round(2πxN​)(1)   由于调用 r o u n d round round取整函数增加了函数调用的开销，这里用其他方式实现 r o u n d round round函数的功能。若 x > 0 x>0 x>0，则 s i g n = 1 sign=1 sign=1;否则， s i g n = − 1 sign=-1 sign=−1。计算 r o u n d ( x N 2 π ) round(\dfrac{xN}{2\pi}) round(2πxN​)等价于计算 ( i n t ) ( x N 2 π + s i g n ∗ 0.5 ) (int) (\dfrac{xN}{2\pi}+sign*0.5) (int)(2πxN​+sign∗0.5)。

  计算 x x x与所对应的最接近其值的正弦表横坐标 x s i x_{s_i} xsi​​的偏差： d = x − s i 2 π N (2) d=x-s_i \dfrac{2\pi}{N} \tag 2 d=x−si​N2π​(2)   由于 c o s ( x ) = s i n ( x + π / 2 ) cos(x)=sin(x+\pi/2) cos(x)=sin(x+π/2)，因而，求 c o s ( x ) cos(x) cos(x)时，实数 x x x所对应的最接近其值的正弦表索引(不一定在 0 至 N − 1 0 至 N-1 0至N−1之间): c i = s i + N / 4 (3) c_i=s_i + N/4\tag 3 ci​=si​+N/4(3)   由于除法运算较慢，上式可以采用移位运算代替： c i = s i + ( N > > 2 ) c_i=s_i +( N>>2) ci​=si​+(N>>2)   将 s i , c i s_i,c_i si​,ci​变换到 0 至 N − 1 0至N-1 0至N−1之间： { s i = s i & ( N − 1 ) c i = c i & ( N − 1 ) (4) \begin{cases} s_i=s_i \& (N-1) \\ c_i=c_i \& (N-1) \\ \tag 4 \end{cases} {si​=si​&(N−1)ci​=ci​&(N−1)​(4)   其中， & \& &表示按位与。   由于 s i n ( x ) = s i n ( x s i + d ) = s i n ( x s i ) c o s ( d ) + c o s ( x s i ) s i n ( d ) sin(x)=sin(x_{s_i}+d)=sin(x_{s_i})cos(d)+cos(x_{s_i})sin(d) sin(x)=sin(xsi​​+d)=sin(xsi​​)cos(d)+cos(xsi​​)sin(d)， d d d很小， s i n ( d ) , c o s ( d ) sin(d),cos(d) sin(d),cos(d)用等价无穷小代替： { s i n ( d ) ∼ d c o s ( d ) ∼ 1 − d 2 / 2 (5) \begin{cases} sin(d) \sim d \\ cos(d)\sim1-d^2/2 \\ \tag 5 \end{cases} {sin(d)∼dcos(d)∼1−d2/2​(5)   得到： s i n ( x ) = s i n e T a b l e [ s i ] + ( s i n e T a b l e [ c i ] − 0.5 ∗ s i n e T a b l e [ s i ] ∗ d ) ∗ d (6) sin(x)=sineTable[s_i] +(sineTable[c_i] - 0.5*sineTable[s_i]*d) *d \tag 6 sin(x)=sineTable[si​]+(sineTable[ci​]−0.5∗sineTable[si​]∗d)∗d(6)   特别的，当需要同时计算 s i n ( x ) , c o s ( x ) sin(x),cos(x) sin(x),cos(x)时，为了进一步减小计算量，避免重复计算 d , s s , c i d,s_s,c_i d,ss​,ci​，可以将计算 s i n ( x ) , c o s ( x ) sin(x),cos(x) sin(x),cos(x)写成一个函数。    s i n ( x ) sin(x) sin(x)在 0 0 0处的泰勒展开为： s i n ( x ) = x − x 3 / 6 + x 5 / 120 − x 7 / 5040 + . . . sin(x)=x-x^3/6 + x^5/120- x^7/5040+... sin(x)=x−x3/6+x5/120−x7/5040+...。    c o s ( x ) cos(x) cos(x)在 0 0 0处的泰勒展开为： c o s ( x ) = 1 − x 2 / 2 + x 4 / 24 − x 6 / 720 + . . . cos(x)=1- x^2/2 + x^4/24- x^6/720+... cos(x)=1−x2/2+x4/24−x6/720+...。   采用式(5)的等价无穷小表示时， s i n ( d ) sin(d) sin(d)的误差限约为 ϵ s = d m a x 3 / 6 \epsilon_{s}=d_{max}^3/6 ϵs​=dmax3​/6, c o s ( d ) cos(d) cos(d)的误差限约为 ϵ c = d m a x 4 / 24 \epsilon_{c}=d_{max}^4/24 ϵc​=dmax4​/24, d m a x = π / ( N − 1 ) d_{max}=\pi/(N-1) dmax​=π/(N−1)。   计算 s i n ( x ) sin(x) sin(x)的误差限： ϵ s i n e ≤ ∣ s i n ( x s i ) ∣ m a x ∗ ϵ c + ∣ c o s ( x s i ) ∣ m a x ∗ ϵ s = ϵ c + ϵ s (7) \epsilon_{sine}\le |sin(x_{s_i})|_{max} * \epsilon_{c} + |cos(x_{s_i})|_{max} * \epsilon_{s}=\epsilon_{c}+\epsilon_{s} \tag 7 ϵsine​≤∣sin(xsi​​)∣max​∗ϵc​+∣cos(xsi​​)∣max​∗ϵs​=ϵc​+ϵs​(7)

  同理可计算 c o s ( x ) cos(x) cos(x)的误差限： ϵ c o s i n e ≤ ϵ c + ϵ s (8) \epsilon_{cosine}\le \epsilon_{c}+\epsilon_{s} \tag 8 ϵcosine​≤ϵc​+ϵs​(8)   

  快速正弦函数、快速余弦函数的c语言实现如下：